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Abstract 

We propose a new approach to the transfer operator based analysis 
of the conformation dynamics of molecules. It is based on a statistical 
independence ansatz for the eigenfunctions of the operator related to a 
partitioning into subsystems. Numerical tests performed on small systems 
show excellent qualitative agreement between mean field and exact model, 
at greatly reduced computational cost. 

1 Introduction 

Conformation transitions of molecules reflect the global spatial/temporal be- 
haviour of the system, and in particular occur at much slower timescales com- 
pared to the elementary frequencies of the system. In a small peptide with a 
dozen atoms, the typical scale difference is already a factor ~ 10^; for folding 
transitions in proteins, it ranges between 10^ and 10^^, placing these events well 
beyond the timescales accessible via direct trajectory simulation. 

The transfer operator approach, introduced into MD by Deuflhard et al. in 
their fundamental paper [5] (cf. also [11, 12, 4, 6, 7, 10, 14]), allows to access long 
time effects through short time simulations, at the expense of needing to simulate 
'ensembles' of initial conditions, or mathematically: to compute the evolution of 
densities on space. In the latter approach, dominant conformations and their 
transition rates can be identified via the leading eigenvalues and eigenfunctions 
of a suitable transfer operator. The catch is that to compute the latter, which 
are functions on phase resp. configuration space, the number of computational 
degrees of freedom of a direct discretization grows exponentially in the number 
of atoms. 
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To overcome this problem, our goal in this paper is to propose a mean field 
method for computing eigenstates of transfer operators, whose relationship to 
the exact eigenvalue problem for the transfer operator is reminiscent of that 
of Hartree-Fock theory to the many-particle Schroedinger equation in quantum 
chemistry, which overcomes an analogous 'curse of dimension' problem. The 
mean field model only has the dimensionality of a typical strongly interacting 
subsystem, but the densities of the different subsystems are nonlinearly coupled. 

In this paper our goal is to 

• derive the mean field model, 

• validate it both by establishing exact properties and comparing to simula- 
tions of the full problem in low- dimensional examples. 

The theoretical and computational results appear to us to be extremely promis- 
ing. Our theoretical results include mass conservation, energy conservation, and 
asymptotic correctness in the limit of weak subsystem coupling (see Section 4). 
The numerical tests we performed on small systems show excellent qualitative 
agreement between mean field and exact model even when the coupling is of 
order one, at greatly reduced computational cost (see Section 6). 

A fuller theoretical explanation of this good performance even for order one 
coupling would be highly desirable, but lies beyond the scope of the present paper. 
One important aspect of our simulations appears to be a careful choice of subsys- 
tems which makes at least the potential part of the Hamiltonian non-interacting. 
Remarkably, such a choice is always possible in chain molecules with a standard 
force field (consisting of nearest-neighbour bond terms, third neighbour angular 
terms, and fourth neighbour torsion terms), by working in inner coordinates (i.e. 
bond lengths, bond angles and torsion angles). In these coordinates, subsystem 
coupling only occurs through momentum exchange. But physical considerations 
as well as previous work support the belief that fine details of momentum transfer 
do not play a decisive role in conformation dynamics: 

- Molecules in solution are subject to relentless perturbations of momenta due to 
collisions with solvent molecules, yet conformation dynamics robustly takes place 
under these conditions. 

- Many of the standard models in conformation dynamics involve randomiza- 
tion of momenta. In the Langevin equation, this happens by addition of white 
noise to the momentum equation; in the approach by Schiitte in [11], one thinks 
of molecular conformations as subsets of configuration space, considers only a 
spatial transfer operator, and draws momenta at random from their statistical 
distributions in each computational evolution step. 

A more detailed analysis as well as applications of the mean field model to 
large systems are currently in progress and will appear elsewhere. 
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2 Hamiltonian dynamics and Liouville equation 



In situations when quantum effects can be neglected and no bond-breaking or 
bond-formation takes place, the dynamics of a molecule with atoms moving 
about in can be described by a Hamiltonian of form 

H{q,p) = \p-M{qr'p + V{q), (1) 

where (g, p) G M^"^, the mass matrix M is a positive dxd Matrix, and V : M'' ^ M 
is a potential describing the atomic interactions. 

In the case when all degrees of freedom are explicitly included and carte- 
sian coordinates are used, we have d = 3N (where is the number of atoms), 
q = (gi, . . . , gAf) G M^^, p = {pi, . . . ,pn), and M = diagimilsy^s), where qi e R^, 
Pi e M^, rrii > are the position, momentum, and mass of the z*^ atom. In this 
paper we work with the more general form (1), in which the kinetic energy is a 
quadratic form of p depending on q. This form arises when inner coordinates are 
used, which will play an important role below. For an A^-atom chain molecule, 
the latter consist of the (A^ — 1) nearest neighbour bondlengths Vij, the (A^ — 2) 
bond angles 9ijk between any three successive atoms, and the (A^ — 3) torsion 
angles (pijki between any four successive atoms. In order to accurately model con- 
formation changes, V will have to contain at least nearest-neighbour bond terms 
Vij{rij), third neighbour angular terms Vijk{Oijk), and fourth neighbour torsion 
terms Vijke{4>ijke)- In practice the potentials could either come from a suitable 
semi-empirical molecular force field model or from ab-initio computations. 

The Hamiltonian dynamics takes the form 

OH 

q= -Q^iq,p) = M{q)~'p, (2a) 

P = -^(^'P) = ■ Miqr'p^ - VV{q). (2b) 

It will be convenient to denote the phase space coordinates hy z = {q,p) G M^'^ 
and the Hamiltonian vector field by 

/ := i'H ], (3) 



V dq y 

so that (2) becomes 

z = fiz). (4) 

The Liouville equation associated to (2) describes the transport of a passive scalar 
u by the flow associated to eq. (2): 

dtu + f- = 0, (5) 
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where u = u{z,t), u : M^'^ x M ^ M. Because the Hamiltonian vector field / is 
divergence-free, eq. (5) can be written equivalently in the form 

dtu + divz{u f) = 0. (6) 

Since eq. (6) preserves both positivity of u and the total integral J^2d u{z,t) dz, 
it defines an evolution on the space of probability densities 



u e lVr2<^) |m > 0, 



u 



Physically it can be interpreted as an evolution equation for "ensembles" of initial 
data under (2). The evolution law (2) for "sharp" initial data can be recovered as 
a special case: z{t) is a solution to (2) if and only if the delta function transported 
along this solution, u{z,t) = Sz(t), is a solution to (6). 

Finally we note that eq. (6) preserves the (expected value of) energy. 



E{t) := j H{z)u{z,t) dz. 



This is because by an integration by parts 

-E{t) = j H{z)[- div{u{z,t)f{z))^ dz = j V H{z) ■ f{z)u{z,t) dz 
and the inner product VH{z) ■ f\z) vanishes for all z, due to (3). 
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3 Molecular conformations and almost invari- 
ant sets 

A conformation of a molecule - as we understand it [5, 11, 12] - is given by an 
almost invariant (or metastahle) subset of configuration space. 

Roughly speaking, an almost invariant set of a discrete dynamical system 
S : X ^ X is a, subset A (Z X such that the invariance ratio 

fi{S~\A)nA) 

''^^ = — ]Ka) — 

is close to 1, cf. [2]. Here fi denotes a suitable probability measure, typically 
Lebesgue measure or some S'-invariant measure. In our case, S will be the time- 
T-map : M^*^ M^'^ of the Hamiltonian system (2) (More generally, instead 
of a deterministic map, one can consider a stochastic process, cf. [2].) 
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Almost invariant sets can be determined via the computation of the eigenfunc- 
tions at (real) eigenvalues close to one of a certain transfer operator associated 
to S [2, 11]. For example, the Frohenius- Perron operator 

PjJi^A) = fi{S~^{A)), A measurable, 

describes the evolution of measures on phase space. This is a linear operator on 
the space A4c of bounded complex valued measures on X. By definition, < 
WnW and thus the spectrum of P is confined to the unit circle. Eigenmeasures 
fi{A) = fL{S~^{A)) at the eigenvalue 1 are invariant measures. If 7^ G A^c is 
an eigenmeasure of P at the eigenvalue A, then 

A/i(X) = Pfi{X) = fi{S-\X)) = /i(X) 

and thus fi{X) = for A 7^ 1. In particular, if A < 1 and fi are real, then 
there are two positive real measures /i"*", /i~ such that fi = fi~^ — fi~ [Hahn- Jordan 
decomposition). If n is normalized such that = /i"*" + /i~ is a probability 
measure then [2] 

Pl^l(A+)+p|^l(yl-) = A + l 

where A"*" = supp^jX^) and A~ = supp{fi~). Consequently, if A < 1 is close to 1 
then both p\^\{A'^) and p|^|(A~) are close to 1 and thus almost invariant. 



Transfer operators. For measures which are absolutely continuous one can 
equivalently consider P on L^(X, C). Since the flow of (2) is a volume- 
preserving diffeomorphism, this operator takes a particularly simple form in our 
case^: 

P^u = uo (7) 

which is the time-T-map of the Liouville equation (5). Note that for an arbitrary 
function (7 : M — >■ [0, 00) of the Hamiltonian, the function u{z) = g{H{z)) satisfies 
Vzu{z) = g' {H {z))V zH (z) . Thus f-VzU = and u, normalized s.t. / u{z) dz = 1, 
is an invariant density. Of particular interest is the canonical density 

h{z) = Cexp{-(3H{z)), (8) 

C = J exp{—f]H{z)) dz, where (3 = l/{kT) and k is Boltzmann's constant. This 
density describes the distribution of a (constant) large number of molecules at 
temperature T and of constant volume. Note that we can write 

h{z) = h{q,p) = Cp(g)exp (^-^p . M-\q)p^ exp (-/3^(g)) =: hp{q,p)h,{q), 

where Cp{q) and Cg are chosen such that J hp{q,p)dp = 1 for each q, and J hg = 1. 
^we write P"^ in order to stress the dependence of the operator on the integration time T 
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Spatial transfer operator. As mentioned, molecular conformations should 
be thought of as almost invariant subsets of configuration space. Schiitte [11] 
introduced a corresponding spatial transfer operator by averaging (7) over the 
momenta: Let h G L^(M^'') be an invariant density of (7) with h[q,p) = h{q, —p), 
let h{q) = J h{q,p) dp and consider the operator 

S^uiq) = j;^J^ (7r,,$-^(g,p)) hiq,p) dp, (9) 

where 'n'q{q,p) = q is the canonical projection onto configuration space. Schiitte 
[11] showed that under suitable conditions, the spatial transfer operator is self- 
adjoint and quasi-compact on an appropriate weighted space. 



Transition probabilities. A key quantity of interest are transition probabil- 
ities from one region of space into another. The transition probability from a 
region Bi C M^*^ into another region Bj C ffi^*^ in phase space is given by the 
volume fraction of those initial data in B^ which end up in Bj at time T, 



m_ m{^-^{B,)nB,) _ 1 

^) ^^'^ 



where m denotes 2(i-dimensional volume. Using the expression based on the 
transfer operator, similarly transition probabilities between subsets of configu- 
ration space can be defined. Typically, these quantities are sought for large T, 
but only short time simulations are numerically feasible. However, p-j can be 
approximated for large T by repeated matrix-vector multiplications, requiring 
short time evaluations of the fiow only, cf. [3] . 



4 Mean field approximation 

The problem with eq. (6) as it stands is that it is amenable to a direct numerical 
treatment only for very small systems, due to the exponential scaling of the 
number of computational degrees of freedom with particle number. If the phase 
space of each atom is approximated by a i^-point grid, such that the solution 
to (6) at time t becomes a vector in M'^, then the corresponding grid of the A^- 
particle system has K'^ gridpoints and the solution at time t becomes a vector in 
. Our proposal to address this problem is partially inspired by Hartree-Fock- 
and density functional theory methods in quantum chemistry, which allow to 
overcome a related complexity problem for the A^-particle Schrodinger equation. 
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Partitioning into subsystems. Starting point is an, for the moment arbi- 
trary, partition of phase space coordinates z = {q, p) into subsystem coordinates: 



TV 



i=l 



where Pi is the vector of momentum coordinates corresponding to the position 
coordinates g^. Let /j = ^|^, "ff")- Then eq. (2) can be re-written as 

Z, = Mz), 2 = l,..,iV. (11) 

Given a phase space density u{zi, . . . , z^, t) which depends on the phase space 
coordinates of all the subsystems, we introduce reduced densities for each sub- 
system i = 1, . . . , A^, as follows: 

Ui{zi,t) := / u{z,t)dzi, (12) 

J^2{d-di) 

where here and below Zi denotes the coordinates {zj)j^i. Note that, by the 
normalization of u, 

Ui{zi,t) dzi = l, 

i 

that is to say the probability that the i*^ particle has some position and momen- 
tum is one. Next we calculate the exact time evolution of the subsystem densities 
Ui. By (6) and (12), 



dtUi{zi,t) = - J dwz{u{z,t) f{z)) dzi = - I div^X'^{z,t) fi{z)) dzi 
= -div^^ / u{z,t) fi{z)dzi. 



Using (12) this can be rewritten in the form of a single-subsystem transport 
equation, 

dtUi{zi,t) = -div;,^(ui{zi,t)fi{zi,t)^, (13) 
with underlying vector field 

m,t)^l^^^^^. (14) 

J u[z, t) dzi 

Note, however, that (13), (14) is not a closed system, since the vector field fi 
depends on the full density u, not just the subsystem densities Uj (j = 1, .., A^). 
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Mean field model. To close the system (13), (14), we now make the statistical 
independence ansatz 

u{zi, . . . ,ZN,t) = Ui{Zi,t) ■ ■ ■ Un{zn, t). (15) 

The vector field fi{zi,t) in (13) then simplifies to 

fi{z„t)= [ f,{z)l\u,{z,,t)dz,=-. ffh]{z,,t). (16) 

jR2(d-d,) L J 

We call the system of eqs. (13), i = 1,. . . ,N, with /j given by (16) the mean 
field approximation to the Liouville equation. Note that it is a system of 
coupled nonlinear partial integrodifferential equations on the lower-dimensional 
subsystem phase spaces M^* , whereas the original Liouville equation was a linear 
partial differential equation on M^'^, d = X^j(2(ii). Physically, each subsystem 
can be pictured, at time t, as experiencing the force of the ensemble of the other 
subsystems in their "typical" states at time t. 

We record some basic properties of the mean field approximation. 

1. The total densities J Ui{zi,t) dzi are conserved. This is immediate from 
the conservation law form dfUi + div {uif^^) = 0. Thus we may assume 
J Ui{zi,t) dzi = 1 for all i and t, and continue to interpret the Ui as proba- 
bility densities. 

2. For non-interacting subsystems, i.e. 

N 



^^^^ = Y^i^p^ ■ M^ii^y'pi + 

i=l 



it is exact, that is to say if the Ui{zi, t) evolve via (13), (16), then the product 
ui{zi, t) ■ ■ -UNizN, t) solves the original Liouville equation (6). This follows 
from the fact that in this case, 

OH OH 

— = Mi{qi)-^pi, -— = —Q^\Pi ■ Mi{qi)-^p, - VV^{qi) 

and so fi{z) depends only on 2j, as a consequence of which the exact vec- 
torfield /j, the reduced vectorfield fi in (14), and the mean field vectorfield 
(RHS of (16)) all coincide, regardless of the mean field ansatz (15). 

For given Uj, j ^ i, the transport equation for Ui has the form of a Liouville 
equation coming from an associated time-dependent subsystem Hamilto- 
nian, 

Hf{q„p„t) = I H{q,p)\{u,{qj,p,,t)dz,, (17) 
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that is to say 



^ Hf{p,,q,,t) 



In particular, is divergence-free. Note that time-dependence of the 
effective subsystem Hamiltonian enters only through time-dependence of 
the Uj, j 7^ i. 



4. The total energy 



E{t) := / H{z)ui{zi,t)---UN{zN,t)dzi 



■ dz 



N 



is conserved. To see this, calculate using the evolution equation for the Ui 
and an integration by parts 
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Hiz) 



N 



J2^izht)Ylue{ze,t) 



Hiz) 



i=l 
N 



dz\ ■ ■ ■ dzis[ 



^ - div^^ {ui{zi, t)f^\zi, t)) JJ Ui{zi, t) 



dz\ ■ ■ ■ dz 



N 



N 

E 

i=l 



V^^H{z)Y[ue{zi^'t) dzi 



■ Ui{zi,t) /f {zi,t) dz, 



But the term in square brackets equals V ZiH^\zi,t), and hence its dot 
product with ff^\zi,t) vanishes for all Zi and t, on account of (18); conse- 
quently iE{t) = 0. 

Property 2. contains useful information regarding how the, up to now ar- 
bitrary, partitioning into subsystems should be chosen in practice. In order 
to maximize agreement with the full Liouville equation (5), the subsystems 
should be only weakly coupled. In the case of an A^-atom chain, this sug- 
gests to work with subsystems defined by inner, not cartesian, coordinates (as 
is done in the simulation of n-butane in Section 6.4 below). Namely, in in- 
ner coordinates, at least the potential energy decouples completely for stan- 
dard potentials containing nearest-neighbour bond terms, third neighbour an- 
gular terms and fourth neighbour torsion terms: V{{rij), {OijiS)ijk, {(f>ijke)ijke) = 

ijki) • 

A deeper, and perhaps surprising, theoretical property of the mean field model 
which goes beyond property 2. concerns weakly coupled subsystems. Consider 
a Hamiltonian of form H{z) = Hq{z) + eHi^itiz), where Hq is a non-interacting 
Hamiltonian of the form given in 2., and e is a coupling constant. It can then 
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be shown that the exact subsystem densities Ui obtained from (5), (12) and the 
mean field densities obtained by solving (13), (16) differ, up to any fixed time 
T, only by O(e^), not the naively expected 0(e). This means that the effect 
of coupling between subsystems is captured correctly to leading order (in the 
coupling constant) by the mean field approximation. For a proof of this fact see 
our companion paper [8]. 



5 The mean field transfer operator 

The mean field approximation to the Liouville equation introduced above gives 
rise in a natural way to a mean field approximation of the transfer operator 

PV = ^-,r), 

where u is the solution to the Liouville equation (5) with initial condition u{z, 0) = 
u^{z). We define the mean field transfer operator as 

P^f(uO) = u(-,T), 

where u{z,t) = {ui{zi,t), . . . ,uiy{zN,t)) is the solution to the mean field ap- 
proximation (13), (16), z = 1, . . . , A^, of the Liouville equation with initial data 
u(^, 0) = u°(z) = {u1{zi), . . . ,u%{zn))- Note that the components of P^^ are 
multilinear, i.e. for fixed Uj, the map 

is linear. 



The mean field spatial transfer operator. Consider an ensemble of molecules 
whose distribution in phase space is given by an invariant density h{q,p). We 
would like to describe distribution changes in the position space only (cf. Sec- 
tion 3). Following [11] we define the spatial transfer operator 



S^w{q) = j {w{q)h{q,p)) dp, (19) 



where h is the conditional density of p for a given q, i.e. 

h{q,p) 



h{q,p) 



Jh{q,p)dp' 



Now we define the spatial transfer operator corresponding to the mean field sys- 
tem. The distribution of the i-th subsystem is given by 



hi{zi) = / h{z) dz. 



10 



The conditional density of the Pi for given is 



h (q- p ) = ^^'^^"^^^ 

J hi{qi,pi)dpi' 

We therefore define the mean field spatial transfer operator as 

5^f(w°) = w(-,T), (20) 
where w°(g) = (ty?(gi), • • ■,w%{qN)), w(g,t) = {wi{qi,t), . . .,WN{qN,t)) and 

.(g.,T) = I K(u°)]^(g„p,)dp. (21) 



with u^[qi,Pi) = Wi{qi)h{qi,pi), i = 1, . . . , N. Again, the components of S'^j are 
multihnear, i.e. for fixed Wj, the map 

Sl,{^^)wr.= KM]. (22) 

is hnear. 



Mean field eigenfunction approximation. Our goal is to compute eigen- 
functions of the (full) spatial transfer operator (19). As approximations, we are 
going to use products of eigenfunctions (at the leading eigenvalues) of the linear 
component maps (22). In computing these, we fix Wj to the invariant density of 
S'j^£. This approach is motivated by the following observation for non-interacting 
subsystems: Let Si, S2 X ^ X he two maps and Pi, P2 : the associ- 

ated Frobenius- Perron operators. For the product map S = Si ^ S2 '■ X"^ X"^, 
S{xi,X2) = {Si{xi) , S2{x2)) , the Frobenius-Perron operator is given by 

{Ph)ixi,X2) = {Plhi){Xi) ■ {P2h2){x2), 

if h{xi,X2) = hi{xi)h2{x2). Now let Pihi = Xihi and P2^2 = ^2/^2 for some 
eigenvalues Ai, A2 G C, then A1A2 is an eigenvalue of P with eigenfunction hih2. 

In order to compute the invariant density, we resort to an iterative procedure 
which is inspired by the famous Roothaan algorithm from quantum chemistry: 
For a given initial guess w° we iteratively compute a sequence w^, k = 0, 1, . . . 
of approximate invariant densities of S'^i by computing the fixed point of each 
linear component map S'j^f(wf), i.e. by computing w'^+^ = {wi~^^, . . . ,w'p~^) such 
that 

SLi^Dw'l^' = w^-'\ t = l,...,N. 
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6 Numerics and examples 



6.1 Ulam's method 

The (spatial) transfer operator is typically considered as an operator on a (suit- 
ably weighted) [p = 1 or 2) space. As such, it is not directly amenable to 
numerical computations. In [13], Ulam proposed the following discretization: 
Let Xn = {Xi, . . . ,Xn} be a disjoint partition of (phase or configuration) space, 
Vn := spanjxi, . . . ,Xn}, where Xi denotes the characteristic function of and 
Qn the projection from onto Vn defined by 

n If 

Qnf ■■=^CiXi with Ci:=—— f dm. 



The discretized (spatial) transfer operator S'^ : Vn Vn is defined to be 

~ QnS ■ 

Note that Sl can be represented by a stochastic matrix, where the matrix entries 
i^n)ij transition rates between the sets Xj and Xj. In other words, 

{S'^)ij is the probability, that 7rg$^(g,p) G Xj, if g G Xj is sampled according to 
a uniform distribution and then p according to h{q, ■). This yields the numerical 
computation of the transition matrix by a Monte-Carlo method: 



k=l 



where the points q^^\ . . . , g'-'^-* are chosen i.i.d. from Xj according to a uniform 
distribution, and the correspondig p^^^ according to h{q^''\-). Observe that the 
Markov-structure of the transfer operator is preserved: and Qn preserve 

integral and positivity, and the Monte-Carlo approximation of the discretized 
transfer operator is a stochastic matrix, too. 

Note that we have to evaluate the flow map several times for each partition 
element Xi, but that these are short time simulations only. The time T merely 
has to be large enough that motion can be observed. In particular, T may be 
much smaller than characteristic times, e.g. for conformational changes; they can 
be several orders of magnitude apart. Thus, the problem of evaluating the flow 
at time T is well conditioned, and we may use low order explicit schemes for the 
time integration. 
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6.2 Complexity 

Let us first investigate the costs of setting up tlie discretized transfer operator for 
an arbitrary subsystem. Using Ulam's method, we need to perform the following 
steps for each partition element Xj: 

• sample g'-'^'^ G Xj and correspondingly p^^\ 

• integrate the mean field system for time T and initial data (p^''\ Q^'^^), 

• project the endpoint onto position space and find the partition element Xi 
it is contained in. 

Using the canonical density for the invariant density h, there is an explicit rep- 
resentation for the momentum distribution h{q, ■) which can be sufficiently well 
approximated by a linear combination of Gaussians. The numerical time integra- 
tion of the initial points requires several evaluations of the mean field vector field 
(18). This in turn requires the numerical evaluation of a 2{d — di) dimensional 
integral. The intergal w.r.t. the pi can be handled analytically and by an apriori 
computation which is independent of the Wj, Pi and g^. Naively, this leaves us with 
a. d — di dimensional integral. However, note that in the case of non-interacting 
subsystems, i.e. fi{z) = fi{zi), fi can be pulled out and the integral reduces to 
1. For systems with small subsystems, i.e. di < d and d small, and in which only 
a fixed and small number of neighboring subsystems interact, the dimensionality 
of the integral is Ylj^^i^j ~ ^(d), where j ^ i means all subsystems j which 
interact with subsystem i. 

The solution of the resulting eigenvalue problems is simple compared with the 
assembling of the discretized mean field transfer operator, particularly since we 
are only interested in the dominant part of the spectrum. Arnoldi type iteration 
methods can be used. 

6.3 Example: a simple 2d-system 

For (gi,g2) £ consider the potential 



with q; = 3, cf. Figure 1. 

Figure 2 shows the eigenf unctions at the leading eigenvalues of the full spatial 
transfer operator. 



V{qi,q2) 
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Figure 1: Potential (23) of the simple 2d-system. 




% q, 



Figure 2: Eigenfunctions at the second, third and fourth largest eigenvalue of the 
full spatial transfer operator (from left to right). 
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-2-1012-2-1012 
q,,resp. q,,resp. ci. 

Figure 3: Invariant density (left) and eigenfunctions at the second eigenvalue of 
the two components of the mean field spatial transfer operator. 

The mean field approximation to the Liouville equation in this case reads 
explicitly 

Figure 3 shows the invariant density w = (1^1,1^2) (after ten iterations of the 
Roothaan type iteration) as well as the eigenfunctions vi and V2 at the second 
eigenvalue of the two linear component maps «S'j^f(wj), z = 1,2. Figure 4 finally 
shows the functions 



Wi ® fi, V2 ® W2 and Vi ® f2 



which serve as approximations to the eigenfunctions from Figure 2. The qual- 
itative structure of the eigenfunctions of the full spatial operator is quite well 
captured. 




Figure 4: Mean field approximations to the eigenfunctions from Fig. 2. 
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6.4 Example: a model of n-butane 

As a more realistic test case we analyse the n-butane molecule CH3-CH2-CH2- 
CH3, cf. Figure 5. More precisely, we consider a united atom model [1] of this 




Figure 5: Cis- and trans- configuration of ra-butane. 

molecule, viewing each CH3 resp. CH2 group as a single particle. Consequently, 
the configuration of the model is described by six degrees of freedom: three bond 
lengths, two bond and one torsion angle. In order to be able to compare the 
results of our mean field approach to a transfer operator based conformational 
analysis on the full configuration space, we further simplify the model be fixing 
the bond lengths at their equilibrium tq = 0.153 nm. For the bond angles we use 
the potential 

V2{e) = -keicosi9-9o)-l) (25) 
with kg = 65^ and Oq = 109.47 ° and for the torsion angle we employ 

1/3(0) = /sT^ (1.116 - 1.462 cos 1.578 cos^ + 0.368 cos^(/) (26) 
+3.156 cosV + 3.788 cos^ (p) , 

with = 8.314^, cf. Figure 6, see also [9]. We fix nip = 1.672 ■ lO-^^g 
as the mass of a proton and correspondingly mi = 14 nip and ni2 = 15 nip 
as the masses of a CH2 and CH3 group, respectively. With q = (^1, 6^2,0)^ G 
[0, tt] X [0, tt] X [0, 27r] denoting the configuration of our model, the motion of our 
system is determined by the Hamiltonian 

H{q,p) = lp''M{q)-'p + V{q), (27) 

with V{q) = V2{qi) + V2(g2) + ^3(^3) and the mass matrix M(g). The latter 
is computed by means of a coordinate transformation q q{q) to cartesian 
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Figure 6: Potential of the torsion angle. 

coordinates q G M^^ for the individual particles, assuming that there is no external 
influence on the molecule and its linear and angular momentum are zero: We have 

q = Dq{q)q 

and consequently 

M{q) = Dq{qyMDq{q), 

where M denotes the (constant, diagonal) mass matrix of the Hamiltonian in 
cartesian coordinates. 

Results for the full operator 

Figure 7 shows the eigenvectors of the full spatial transfer operator for system 
(27) at the two largest eigenvalues ^ 1 (computed on a 32 x 32 x 32 grid using 
32 sample points, T = 0.5 ■ 10~^^ s integration time, realized by 10 steps of the 
explicit Euler method with step size T/10). 

Results for the mean field approach 

We decompose the model into three subsystems, i.e. each configuration vari- 
able is treated separately. The Roothaan iteration is initialized with iff(gi) : = 
C.e-/3^>fe), i = 1,2,3, where p is the inverse temperature corresponding to 300 
K and Ci is a corresponding normalizing factor. Figure 8 shows the mean field 
approximations to the two eigenvectors in Fig. 7, namely the products 

® v^,2 and we^ ® w</,,3, 

where is the ^2-factor of the invariant density of the mean field system and 
^</),2, w</>.3 are the eigenvectors at the second and third largest eigenvalue of the 
(j) subsystem, respectively (after ten iterations of the Roothan type iteration). 
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Figure 7: Eigenvectors of the full spatial transfer operator at A2 = 0.985 
and A3 = 0.982 (right). Shown is a slice at qi = 61 = 7r/2. 




Figure 8: Mean field approximations to the two eigenvectors shown in Fig. 
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Clearly, the qualitative structure of the eigenvectors of the full operator (cf. Fig. 7) 
is well captured by the mean field approximation. 

7 Conclusion and Outlook 

We introduced a new approach to the transfer operator based conformational 
analysis of molecules in this paper. The central idea is a mean field approach 
based on a statistical independence ansatz for the eigenfunctions of the opera- 
tor. The principal motivation for such an ansatz lies in its linear, as opposed to 
exponential, scaling of the number of computational degrees of freedom with the 
number of atoms, provided the subsystem size stays fixed. We established basic 
theoretical properties including mass and energy conservation and asymptotic 
correctness in the limit of weak subsystem coupling. In numerical tests on small 
systems, the mean field model is seen to provide a remarkably accurate represen- 
tation of the true eigenfunctions. Applications to larger systems are currently in 
progress and will be discussed elsewhere. 
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